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Abstract 



We offer a new proposal for the Monte Carfo treatment of many-fermion 
systems in continuous space. It is based upon Diffusion Monte Carlo with 
significant modifications: correlated pairs of random walkers that carry oppo- 
site signs; different functions "guide" walkers of different signs; the Gaussians 
used for members of a pair are correlated; walkers can cancel so as to con- 
serve their expected future contributions. We report results for free-fermion 
systems and a fermion fluid with 14 ^He atoms, where it proves stable and 
correct. Its computational complexity grows with particle number, but slowly 
enough to make interesting physics within reach of contemporary computers. 
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Monte Carlo methods have provided powerful numerical tools for quantum many-body 
physics [|1],|2[. They include methods such as Green's function Monte Carlo (GFMC) 0, 
Diffusion Monte Carlo (DMC) [g, or Path Integral Monte Carlo (PIMC) that are ca- 
pable of giving, at least for moderate size bosonic systems, answers with no uncontrolled 
approximations. Such accurate treatment of fermionic systems has been made vastly more 
difficult by a "sign problem." Progress in the application of Quantum Monte Carlo methods 
to condensed matter physics, to electonic structure, and to nuclear structure physics has 
been impeded for years by the lack of exact and efficient methods for dealing with fermions. 

This paper offers a new proposal for solving many-fermion systems by an extension of 
DMC. In the systems we have studied, the signal-to-noise ratio of the Monte Carlo esti- 
mates are constant at long imaginary times, by contrast to the behavior of ordinary DMC 
where they decay exponentially 0. Except for the use of a short-time Green's function, 
no approximations- physical, mathematical, or numerical- are made. The effect of a finite 
interval of imaginary time is easily controlled, or may be eliminated entirely. 

It is no surprise that Monte Carlo methods can solve the many-body Schrodinger equation 
in imaginary time for bosonic systems. Let R denote all coordinates of an A^-body system, 
and V{R) be the potential at R. 



This equation also describes the diffusion of an object (a "random walker") in a 3A^ di- 
mension space in which the potential V{R) serves as a generalized absorption rate. Because 
the potential in physical problems can be unbounded from above and below, a direct simula- 
tion of that diffusion, although straightforward, will be inefficient. Some form of importance 



this is a technical device for accelerating the convergence; in our new method it becomes an 
essential feature. 



construct a random walk. A simple version is as follows: Using a fixed step in imaginary time. 




(1) 



sampling transformation has been found to be highly useful. In the standard DMC 



DMC uses an "importance" or "guiding" function iPq{R) and a trial eigenvalue Et to 
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5r, a walker at i? is (a) moved to R+ StV \n ipGiR)] (b) then each coordinate is incremented 
by an element of a random vector U, a Gaussian with mean zero and variance 6t; finally, 
(c) each walker is turned into M walkers with < M >= exp |5r Ex — HiI)g{R) / i^ciR) 
where H is the Hamiltonian. 

The resulting random walk has expected density 

/(i?,r) = V^G(i?)E«A.exp[(ET-^fc)r]0fe(i?) (2) 

k 

where 4>k{R) are eigenf unctions of H with eigenvalues E^, and are expansion coefficients. 
The limit of f{R, t) for large r is dominated by the eigenfunction (pQ having the lowest 
eigenvalue Eq. 

We alter the structure of DMC in the following ways: (1) In order to represent an anti- 
symmetric wave function that is both positive and negative, we introduce walkers, {R^, Rm}^ 
that respectively add or subtract their contributions to statistical expectations. The com- 
putation now involves ensembles of pairs of walkers carrying opposite signs. (2) Two distinct 
functions, '4^g{R) are used to guide the ± walkers. (3) The Gaussians for the paired 
walkers are not independent; rather U~ is obtained by refiecting in the perpendicular 
bisector of the vector i?+ — R^ . Finally (4) the overlapping distributions that determine 
the next values of i?^ are added algebraically so as to allow positive and negative walkers 
to cancel, while preserving the correct expected values. 

In a Monte Carlo calculation of this kind, we "project" quantities of interest by calculat- 
ing integrals weighted with some trial function, say tpriR)- In DMC the energy eigenvalue, 
Eq, can be determined from: 

° fMR)MR)dR ^ MRm) 

replacing integrals by sums over positions of the random walk. 
In our modified dynamics, Eq.(H) now takes the form 
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If {-Rm} and {R^} follow the same dynamics using the same V'g, then the expected values 
of numerator and denominator of Eq. (H) decay exponentially at large r. Some correlation 
among walkers is essential. This observation is reinforced by noting that the Pauli principle, 
which demands that fermion wave functions be antisymmetric, is a global condition, and 
cannot be satisfied by independent walkers. Put another way, it will be necessary to have 
dynamics that distinguish between walkers that carry different signs. These motivations 
underlie aspects (2) and (3) of our method outlined above. 

A second aspect of the difficulty in treating fermion systems is that the density that one 
obtains naturally from a random walk is the symmetric ground state. In order for Eq. 
to have an asymptotically bounded signal-to-noise ratio, walkers of opposite signs must be 
able efficiently to cancel. This underlies the need for modification (4) given above. The 
need for some degree of cancellation has been a theme of previous research starting with the 
work of Arnow et al. 0. The need for distinct dynamics for positive and negative walkers 
was stressed in p|. That these two aims could be accomplished by appropriate correlation 
among walkers was pointed out by Liu, Zhang, and Kalos 0. The use of distinct guiding 
functions is new and serves as the connection among the different algorithmic ideas that 
enables the treatment of general potentials. 

Stable results can be obtained using correlated pairs only. Correct results are ensured 
when the dynamics have the property that the random walk for either member of the pair is 
the same as that of a single free walker, except when they cancel, a condition satisfied here. 
The expectations of Eqs.(^) and (^) are linear in the walker densities and are unchanged 
by correlations. Furthermore, we have devised a method of canceling opposite walkers that 
also preserves these expectations. 

— # — * 

Let ^pAiR) be a trial function for the fermionic state. Let fsiR) be some approximation 
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to the symmetric ground state wave function of the same Hamiltonian. Define: 



^P^iR) = ^cfliR) + c^^\{R) ± cvj,{R) (5) 

The following properties of these two functions are significant: (a) they are positive; (b) 
when c is small, they are dominated by ^psi so that opposite walkers behave similarly; (c) 
ipQ transforms under an odd permutation V as follows: 

^UvR) = rG{R)- (6) 

As mentioned above we modify simple DMC in several ways. The "drift" is applied in the 
usual way to walkers assumed to be at i?^, using the two guiding functions: 

R+ = R+ + 6TV\n'^^{R+) 

(7) 

R- = R- + STV\n^aiR~)- 

— * 

Diffusion of the walkers, however, is carried out in a correlated way: let be a vector of 
3A^ Gaussian random variables each of mean zero and variance 6t. New trial positions R^ 
are now given by 

R+ = R+ + u+. R- = R-^ ij-^ (8) 

where the random vector U~ is obtained by reflection in the perpendicular bisector of the 
vector _R+ — R^ as described in (4) above. 

Walker densities can be subtracted by computing their chances of arrival at a common 
point, but because they have different guiding functions, they do not exactly cancel. The 
analysis of "forward walking" |I|,|6|,p!0| allows one to determine the expected future contribu- 
tion of a walker to any projected quantity. Thus, we can compute the change in expectations 
when a pair meets. For this change to be zero [O , a positive walker at R^ must survive to 
the next time step with probability 

P+{R+-R+,R~) = 

(9) 

max 0, 1 



B-{Rt^ 


R )G(^R^ 


-R-)^MRt) 


BHR^ 


R^)G[R^ 


~R^)iJaiR^) 



where 



G{R' - R) 



exp[-(#-i?)V(25r)] 



(10) 



(27r5r)3^/2 

is the Gaussian density used in DMC The branching factors, B^{R\R^) and B^{R\R^) are 

H^p^iR) 



B^{R) = exp <^ 6t 



T 



(11) 



An analogous expression is used for negative walkers. 

An isolated walker may appear as a result of different branching factors at {R^} and 
if, with probability one half, one generates a walker of opposite sign by interchanging 
the coordinates of two like-spin particles, then a pair is reconstituted that preserves future 
expectations. 

To determine the energy, we use the estimator of Eq. (p. A sharp indication of the 
stability of the calculation is the behavior of its denominator 



V 



(12) 



>+(i?+) 

In a naive calculation, P decays to zero in an imaginary time of order Tc = 1/ {Ea~Es) where 
Ea and Es are the fermion and boson energies. A stable method will show T> asymptotically 
constant. 

Although a system of free fermions in a periodic box is analytically trivial, it presents 
an exigent test of this method. For this system, the lowest symmetric state is constant, and 
the exact fermionic wave function is a determinant of plane waves. We use p = 0.5 and set 



^^{R) = ^J1 + c^vI{R)±cva{R) 
where ipA is a Slater determinant of one body orbitals Xf", of the following form: 



(13) 



An 



exp 



ik- \ri + XbY,v{ 



(14) 



The parameter \b controls the departure of the nodal structure of this function from the 
exact shape. The fact that these functions are modulated only a little from a constant by 
if A means that the polarization of the population of plus and minus walkers is small. 



In table I we report the results obtained for periodic systems of 7, 19, and 27 free 
fermions. The results agree with the analytic eigenvalues within the Monte Carlo estimates 
of the standard error. It has been conjectured that the computational complexity of Fermion 
Monte Carlo calculations will grow as A^!, where N is the number of particles in the system. 
Since (27!/7!) = 2.16 xlO^'', a calculation with 27 or even 19 bodies would be impossible 
were that conjecture to be true. 

We have also applied this algorithm to a system of 14 ^He atoms in a periodic box at 
equilibrium density, p = 0.0216A~^. Energies are expressed in Kelvins, and lengths in A. 

With interatomic potentials that have a hard core, we may use the same function ipA as 
for free fermions, but also need a Jastrow product. With 



In Fig. 1 we plot the cumulative denominator as a function of imaginary time for a 
typical run. As can be seen, the fundamental stability of the method is well demonstrated. 
Fig. 2 shows the decay of the same denominator, when the method is made unstable by 
setting c = 0. 

Table II exhibits the eigenvalues of various runs with our method applied to the periodic 
system with 14 ^He atoms. They are all consistent, and yield a weighted average of - 
2.2558(39). The run marked (b) is a continuation of the run labeled (a) separated by a 
long run with a longer time step. As a whole, including such continuations, the longest 
aggregate sequence comprises a total imaginary time of 1830 K~^. Using a total system 
energy difference of 20 K (as we have measured), that corresponds to 3.6 x 10^ fermion 
decay times. An alternative measure of the length of the run, suggested by David Ceperley 
|T^, is the ratio of the rms diffusion length of a particle to the mean spacing between 
particles. For this sequence of runs, that ratio is 19. Thus the observation of stable values 
of the sums in Eq. (^ is significant. 



= ¥^s{R) = l[exp[-{b/rij)% 



(15) 



i<j 



the guiding functions now have the form: 




(16) 
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Space limitations preclude a complete description of the other checks that we have made 
that the results for ^He are correct: they include a fixed node calculation of exactly the same 
model problem, which yielded an eigenvalue of -2.08(1) K. A transient estimate (cf. Fig. 
3), relaxing from the fixed node, is consistent with our result (shown as the dashed line.) 
Analysis of the results in Fig. 2 leads to a fermion-boson energy difference of 1.434(35) K 
per particle. This agrees well with a direct calculation of the energy of a 14-body mass-3 
boson system that gave -3.68(1) K. 

By construction, the method proposed here introduces no approximations other than 
that of the short imaginary-time Green's function. In other words, if the results are stable, 
then they are correct. Although we have have not yet proved the stabihty of the method 
(i.e. that the long-term average of the denominator of the eigenvalue quotient is not zero), 
we believe that we have convincingly demonstrated the stability. Perhaps the most impor- 
tant conclusion that we may draw is that the "sign problem" of Fermion Monte Carlo for 
continuous systems is not intractable; the search for elegant computational methods in this 
and related applications is justified. 
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TABLES 



N 




E 




■'-'ex 


7 




2.91285(49) 




2.912712 


19 




2.760(25) 




2.757454 


27 




2.796(30) 




2.763316 


TABLE I. 


Energies and errors for 


a periodic system of N free fermions. The 


analytic result is 












b 


c 




As 


E 


1.15 


0.025 







-2.251(08)a 


1.145 


0.025 







-2.258(17) 


1.135 


0.025 







-2.257(10) 


1.15 


0.016 







-2.246(20) 


1.15 


0.010 







-2.250(19) 


1.15 


0.025 







-2.2559(84)b 


1.15 


0.025 




-0.05 


-2.249(12) 


1.15 


0.025 




0.05 


-2.268(10) 



TABLE II. Energies and errors for a periodic system of 14 ^He atoms 
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FIGURES 

FIG. 1. Cumulative denominator of energy quotient: a stable calculation 

FIG. 2. Denominator of energy quotient: an unstable calculation. The smooth curve is 
exp(— 20.08r), a fit to the data. 

FIG. 3. Relaxation of eigenvalue from fixed node: a transient calculation. 
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http://arXiv.org/ps/cond-mat/0001193vl 



This figure "Mel2.jpg" is available in "jpg" format from: 
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